LiDAR Tutorial using R

Introduction

Packages

library(lidR)# for working with  LAZ files
                            library(sf) # for working spatial class
                            library(raster)# for working with raster
                            library(rayshader) # for 3D viz
                            library(rgl) # for interactive plots

Area of Interest

motte_point = st_read("data/canmore_point.shp", quiet = TRUE)
                            
                            motte_buffer = st_buffer(motte_point,dist = 50)

LAZ

Next, I am reading the LAZ files and clipping the point cloud to the extent of immidiate area around the motte.

#read laz files
                            las = readLAS("data/NX6055_4PPM_LAS_PHASE3.laz")
                            
                            motte = clip_roi(las, motte_buffer)
                            #plot lidar point cloud
                            plot(motte, bg = "white")
                            
                            rglwidget()

grab and rotate the model !

Digital Elevation Model

In this step I am running standard algorithms from LiDR package to compute Digital Surface Model (DSM) and Digitial Terrain Model (DTM).

rgl.clear()
                            # create dsm and dtm
                            dsm = grid_canopy(motte, 0.5, pitfree())
                            
                            # assign coordinate system
                            crs(dsm) = CRS("+init=epsg:27700")
                            
                            # create dtm
                            dtm = grid_terrain(motte, 0.5, algorithm = knnidw(k = 6L, p = 2))
                            
                            # addign coordinate system
                            crs(dtm) = CRS("+init=epsg:27700")
                            
                            par(mfrow = c(1,2))
                            plot(dtm, main = "DTM", col = height.colors(50))
                            plot(dsm, main = "DSM",col = height.colors(50))

Hillshade

# dsm
                            slope_dsm = terrain(dsm, opt = 'slope')
                            aspect_dsm = terrain(dsm, opt = 'aspect')
                            hill_dsm = hillShade(slope_dsm, aspect_dsm, angle = 40, direction = 270)
                            
                            # dtm
                            slope_dtm = terrain(dtm, opt = 'slope')
                            aspect_dtm = terrain(dtm, opt = 'aspect')
                            hill_dtm = hillShade(slope_dtm, aspect_dtm, angle = 5, direction = 270)
                            
                            #plot
                            par(mfrow = c(1,2))
                            plot(hill_dtm, main = "DTM Hilshade", col = grey.colors(100, start = 0, end = 1), 
                                 legend = FALSE)
                            plot(hill_dsm, main = "DSM Hillshade", col = grey.colors(100, start = 0, end = 1))

3D Model

#And convert it to a matrix:
                            elmat = raster_to_matrix(dtm)
                            
                            elmat %>%
                              sphere_shade(texture = "imhof1") %>%
                              add_shadow(ray_shade(elmat, zscale = 0.5, sunaltitude = 30,lambert = TRUE),
                                         max_darken = 1) %>%
                              add_shadow(lamb_shade(elmat,zscale = 0.5,sunaltitude = 30), max_darken = 0.2) %>%
                              add_shadow(ambient_shade(elmat, zscale = 0.5), max_darken = 0.2) %>%
                              plot_3d(elmat, zscale = 0.5,windowsize = c(1000,1000),
                                      phi = 40, theta = 180, zoom = 0.8, fov = 1)
                            rglwidget()

grab and rotate the model !

Another example …

Tantallon Castle And here is the code

library(rayshader)
                            library(raster)
                            library(magick)
                            library(gifski)
                            
                            dem = raster("data/tantallon.tif")
                            
                            n_frames <- 41
                            
                            elmat = raster_to_matrix(dem)
                            
                            zscale <- 0.5
                            
                            waterdepthvalues <- c(0:20, seq(19,0,-1)) / 2
                            
                            phi_values = seq(20,60)
                            
                            waterdepthvalues * zscale
                            
                            length(waterdepthvalues)
                            
                            img_frames <- paste0("drain", seq_len(n_frames), ".png")
                            
                            for (i in seq_len(n_frames)) {
                            elmat %>%
                              sphere_shade(texture = "imhof1") %>%
                              add_shadow(ray_shade(elmat, zscale = 0.5, sunaltitude = 30,lambert = TRUE),
                                         max_darken = 1) %>%
                              add_shadow(lamb_shade(elmat,zscale = 0.5,sunaltitude = 30), max_darken = 0.2) %>%
                              add_shadow(ambient_shade(elmat, zscale = 0.5), max_darken = 0.2) %>%
                              plot_3d(elmat, zscale = 0.5,windowsize = c(1000,1000),
                                      water = TRUE, watercolor = "imhof3", wateralpha = 0.5, 
                                      waterlinecolor = "#ffffff", waterlinealpha = 0.5,
                                      waterdepth = waterdepthvalues[i],
                                      phi = 50, theta = 45, zoom = 0.8, fov = 1)
                              render_snapshot(filename = img_frames[i],
                                              title_text = "Tantallon Castle", 
                                              title_font = "Helvetica",
                                              vignette = TRUE,
                                              title_size = 50,
                                              title_color = "black")
                              rgl::clear3d()
                            }
                            
                            magick::image_write_gif(magick::image_read(img_frames), 
                                                    path = "tantallon.gif", 
                                                    delay = 6/n_frames)
sessionInfo()
## R version 4.0.2 (2020-06-22)
                            ## Platform: x86_64-apple-darwin17.0 (64-bit)
                            ## Running under: macOS Mojave 10.14.6
                            ## 
                            ## Matrix products: default
                            ## BLAS:   /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
                            ## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
                            ## 
                            ## locale:
                            ## [1] en_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/C/en_GB.UTF-8/en_GB.UTF-8
                            ## 
                            ## attached base packages:
                            ## [1] stats     graphics  grDevices utils     datasets  methods   base     
                            ## 
                            ## other attached packages:
                            ## [1] rgl_0.100.54     rayshader_0.15.1 sf_1.0-2         lidR_3.1.4      
                            ## [5] raster_3.1-5     sp_1.4-5        
                            ## 
                            ## loaded via a namespace (and not attached):
                            ##  [1] Rcpp_1.0.7              lattice_0.20-41         png_0.1-7              
                            ##  [4] prettyunits_1.1.1       class_7.3-17            digest_0.6.27          
                            ##  [7] foreach_1.5.0           utf8_1.2.2              mime_0.11              
                            ## [10] R6_2.5.1                evaluate_0.14           e1071_1.7-4            
                            ## [13] highr_0.9               pillar_1.6.2            rlang_0.4.11           
                            ## [16] progress_1.2.2          lazyeval_0.2.2          data.table_1.14.0      
                            ## [19] miniUI_0.1.1.1          jquerylib_0.1.4         magick_2.3             
                            ## [22] rmarkdown_2.10.1        rgdal_1.5-18            webshot_0.5.2          
                            ## [25] stringr_1.4.0.9000      htmlwidgets_1.5.1       munsell_0.5.0          
                            ## [28] shiny_1.4.0.2           compiler_4.0.2          httpuv_1.5.4           
                            ## [31] xfun_0.25               pkgconfig_2.0.3         rgeos_0.5-2            
                            ## [34] htmltools_0.5.1.1       tidyselect_1.1.0        tibble_3.1.4           
                            ## [37] codetools_0.2-16        fansi_0.5.0             crayon_1.4.1           
                            ## [40] dplyr_1.0.2             later_1.1.0.1           grid_4.0.2             
                            ## [43] jsonlite_1.7.2          xtable_1.8-4            lifecycle_1.0.0        
                            ## [46] DBI_1.1.0               magrittr_2.0.1          scales_1.1.1           
                            ## [49] units_0.6-7             KernSmooth_2.23-17      stringi_1.7.3          
                            ## [52] promises_1.1.1          doParallel_1.0.15       bslib_0.2.5.1          
                            ## [55] ellipsis_0.3.2          generics_0.1.0          vctrs_0.3.8            
                            ## [58] iterators_1.0.12        tools_4.0.2             manipulateWidget_0.10.1
                            ## [61] glue_1.4.2              purrr_0.3.4             hms_0.5.3              
                            ## [64] crosstalk_1.1.0.1       parallel_4.0.2          fastmap_1.0.1          
                            ## [67] yaml_2.2.1              colorspace_2.0-2        rlas_1.3.5             
                            ## [70] classInt_0.4-3          knitr_1.33              sass_0.4.0